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Introduction 

The aim of this paper is to propose a practical numerical method to approximate 
the survival function and the moments of the exit time for a piecewise-deterministic 
Markov process thanks to the quantization of a discrete-time Markov chain naturally 
embedded within the continuous-time process. 

Piecewise-deterministic Markov processes (PDMP's) have been introduced by M.H.A. 
Davis in |5] as a general class of stochastic models. PDMP's are a family of Markov 
processes involving deterministic motion punctuated by random jumps. The motion 
depends on three local characteristics namely the flow (f>, the jump rate A and the 
transition measure Q, which specifies the post-jump location. Starting from the point 
X, the motion of the process follows the flow '^{x, t) until the flrst jump time Ti, which 
occurs either spontaneously in a Poisson-like fashion with rate A(<i>(a;,i)) or when the 
flow <i>(x, t) hits the boundary of the state space. In either case, the location of the 
process at the jump time Ti, is selected by the transition measure Ti)) and 

the motion restarts from this new point X{Ti) denoted by Zi. We deflne similarly the 
time 5*2 until the next jump, the next jump time is T2 = Ti -\- S2, the next post-jump 
location Z2 = X{T2) and so on. Thus, associated to the PDMP we have discrete-time 
Markov chains {Zn, T'n)neN5 given by the post-jump locations and the jump times, and 
gN, given by the post-jump locations and the inter-jump times. A suitable 
choice of the state space and the local characteristics <i>, A and Q provides stochastic 
models covering a great number of problems of operations research for example see |1] , 
[S] and the corrosion model presented in this paper. 

Numerical computation of the moments of the exit time for a Markov process has 
been studied by K. Helmes, S. Rohl and R.H. Stockbridge in [5]. Starting from an 
assumption related to the generator of the process, they derive a system of linear 
equations satisfled by the moments. In addition to these equations, they include flnitely 
many Hausdorff moment conditions that are also linear constraints. This optimization 
problem is a standard linear programming problem for which many efficient softwares 
are available. J.-B. Lasserre and T. Prieto-Rumeau introduced in |10j a similar method 
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but they improved the efficiency of the algorithm by replacing the Hausdorff moment 
conditions with semidefinite positivity constraints of some moment matrices. Never- 
theless, their approach cannot be applied to PDMP's because the assumption related 
to the generator of the process is generally not satisfied. In [5] section 33, M.H.A. 
Davis gives an iterative method to compute the mean exit time for a PDMP but his 
approach involves solving a large set of ODEs whose forms are very problem specific, 
depending on the behaviour of the process at the boundary of the state space. Besides, 
and in the context of applications to reliability, it seems important to study also the 
distribution of the exit time. 

There exists an extensive literature on quantization methods for random variables 
and processes. The interested reader may for instance consult [7], [TT] and the ref- 
erences within. Quantization methods have been developed recently in numerical 
probability or optimal stochastic control with applications in finance (see e.g. [T], 
[2], [3] and [llj'l. The quantization of a Markov chain {Qn)neN consists in finding, for 
each n, an optimally designed discretization of the state space of 0„ providing the best 
possible approximation by a random variable 0„ taking its values in a grid r„ of 
finite and fixed size as well as transition measure of the quantized chain (0n)neN- As 
explained for instance in [11 , section 3, provided that the Markov kernel is Lipschitz, 
bounds for the rate of convergence of the quantized process towards the original 
process are obtained. 

In the present work, we consider a PDMP {Xt)f.^f^ with state space E and we present 
approximation methods to compute the moments and the survival function of the exit 
time from a set denoted U C E given the fact that it happens before the iV"^ jump 
of the PDMP denoted by T/y. Roughly speaking, we estimate the moments and the 
survival function for r A Tjy. In our approach, the first step consists in expressing the 
j-th moment (respectively the survival function) as the last term of some sequence 
{Pk,j)k<N (respectively {pk)k<N) satisfying a recursion Pk+i,j = i>{Pk,j) (respectively 
Pk+i = ''PiPk)) specifically built within our paper. 

In this context, a natural way to deal with these problems is to follow the idea 
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developed in [6 namely to write the recursions in terms of an underlying discrete- 
time Markov chain and to replace it by its quantized approximation. The definitions 
of {pk.j)k and (pfe)fe involve some discontinuities related to indicator functions but 
as in we show that they happen with small enough probability. However, an 
important feature that distinguishes the present work from [6 and which prevents 
a straightforward application of the ideas developed within, is that an additional 
important difficulty appears in the definition of the sequences {pk,j)k and {pk)k- Indeed, 
the mapping ip such that Pk+i,j — i^iPkj) andpfe+i = ip{pk) is not Lipschitz continuous. 
One of the main results of this paper is to overcome this difficulty by deriving new and 
important properties of the Markov chain (Z„, r„)„gN, combined to a sharp feature of 
the quantization algorithm. We are able to prove the convergence of the approximation 
scheme. Moreover, in the case of the moments, we even obtain bounds for the rate 
of convergence. It is important to stress that these assumptions are quite reasonable 
with regards to the applications. 

An important advantage of our method is that it is flexible. Indeed, as pointed out 
in U, a quantization based method is "obstacle free" which means, in our case, that it 
produces, once and for all, a discretization of the process independently of the set U . 
Consequently, the approximation schemes for both the moments and the distribution 
of the exit time are flexible w.r.t. to U. Indeed, if we are interested in the exit time 
from a new set C/', it will be possible, provided U' satisfies the same assumptions as 
U, to obtain in a very simple way the moments and the distribution of this new exit 
time. Indeed, the quantization grids are only computed once, stored off-line and may 
therefore serve many purposes. 

The paper is organized as follows. We first recall the definition of a PDMP and state 
our assumptions. In Section [2] we introduce the moments and the distribution problems 
and present recursive methods to solve them. Section |3] contains the main contribution 
of this paper namely the approximation schemes, the proofs of convergence and bounds 
for the rates of convergence. Eventually, two numerical examples are developed in 
Section [4] and the advantages of our approach are discussed in Section [5] 
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1. Definitions and assumptions 



For any metric space X, we denote B{X) its Borel a-ficld and B{X) the set of 
real- valued, bounded and measurable functions defined on X. For a,b G M, denote 
a Ab = min(a, b) and a V 6 = max(a, 6). 

Definition of a PDMP 

In this first section, let us define a piecewise-deterministic Markov process and 
introduce some general assumptions. Let M be a finite set called the set of the modes 
that will represent the different regimes of evolution of the PDMP (M is supposed to 
be a finite space although it could be countable), for each m G M, the process evolves 
in Erm open subset of M''(") (where d: M ^W). Let 



This is the state space of the process {Xt)ti£R+ = {'mt,^t)teR+- Let dE be its boundary 
and E its closure and for any subset Y of E, V denotes its complement. 

Define on E the following distance, for x = (m, ^) and x' = (m',^') € E, 



I otherwise. 

Moreover, for any x G E and Y C E, denote d{x,Y) the distance between the point x 
and the set Y i.e. d{x, Y) = infy^y — y\- 

A PDMP is defined by its local characteristics ($m, ^m, Qm)meM '■ 

• For each m e M, : M''(™) x M ^ K^fC™) is a continuous function called the 
flow in mode m. For all f e R, is an homeomorphism and t ^m{-,t) 

is a group i.e. for all ^ e ]R<^(™), ^m{^,t + s) = $„($„(^,s),i). For all x = 
(m, ^) e E, define now the deterministic exit time from E : 



E = {{m,^),m G M,^ G Em} ■ 




(1) 



e{x) 



inf{t > such that $r»(^,t) G dEm} 



We use here and throughout the whole paper the convention inf = +oo. 
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• For all m e Af , the jump rate : M+ is measurable and satisfies : 

V(m,f) e i?, 3e > such that / \jn{<^yri{i,t))dt < +00. 

Jo 

• For all m G M, Qm is a Markov kernel on {B{E), Em) which satisfies : 

y^eEm, Q,„(i?\{(m,0},0 = l- 

From these characteristics, it can be shown (see [5]) that there exists a filtered prob- 
ability space {fl,J-',J-'t,{Px)xeE) on which a process {Xt)t£M+ is defined. Its motion, 
starting from a point x G E, may be constructed as follows. Let Ti be a nonnegative 
random variable with survival function : 

g-A(x,t) ifo<t<r(a;), 
iit>t*{x), 
where for x = (m,^) e E and t e [0,t*{x)], 

A{x,t) = [ A™($™(f,s))ds. 
Jo 

One then chooses an _B- valued random variable Zi with distribution Qmi', ^miS,, Ti)). 
The trajectory of Xf for i < Ti is : 



PATi > t) 



Xt 



(m,$™(e,t)) ift<Ti, 
Zi iit^Ti. 

Starting from the point Xti — Z\ , one then selects in a similar way Si — T^— T\ the 
time between T\ and the next jump, Z2 the next post-jump location and so on. 
M.H.A. Davis shows (see [5]) that the process so defined is a strong Markov process 
{Xt)t>o with jump times (r„)„gN (with To 0). The process (6„)„eN = (^n,r„)„gN 
where Z„ = Xt^ is the post-jump location and r„ is the rt-th jump time is clearly 
a discrete-time Markov chain. Besides, we denote Sn — Tn — Tn-i and 6*0 = the 
inter-jump times. 

The following assumption about the jump-times is standard (see for example ^ , section 
24) : 

Assumption 1.1. For all {x,t) € E x R+ , E^, ^{n<t}] < 
It implies that Tk — >■ +00 a.s. when k — > +00. 
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For notational convenience, any function h defined on E will be identified with its 
component functions hm defined on Em- Thus, one may write 

h{x) = hmiO when x = (m, ^) € E. 

We also define a generalized flow ^ : E x M+ — >■ E such that 

^{x,t) = (m, when x = (to, ^) G E. 

Notation 

For any function w in B{E), introduce the following notation 

Qw{x) = / w{y)Q{dy,x), Cyj = sup Iwix)] 
Je x&E 

and for any Lipschitz continuous function w in B{E), denote [w] its Lipschitz constant; 

\w{x) ~ w{y)\ 



[w] = sup_ _ 



with the convention — = 0. 



Remark 1.2. For w € B(E) and from the definition of ttie distance on E, one iias 

[w] = SUPmeM[Wm]- 

2. Exit time 

For all m G M, let Um be a Borel subset of Em, let [/ = {(to,^),to G M, ^ G Um}- 
We are interested in the exit time from U denoted by r ; 

T = inf{s > such that ^ C/}. (2) 

Denote fi the distribution of the initial state of the process Zq. Since the present paper 
concerns numerical computations, the following assumption appears natural. 

Assumption 2.1. The process starts in U and eventually leaves it almost surely i.e. 
the support of n is included in U and (r < +oo) = 1. 

The aim of this paper is to provide approximation schemes for its survival function 
and its moments. Our method has a high practical interest because it will provide 
numerical approximations as soon as the process can be simulated. Our approach is 
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based on a recursive computation using the underlying discrete-time Markov chain 
(Z„, T„)„gN. Therefore, we will study r A Tjv rather than r for some iV G N called the 
computation horizon. Indeed, thanks to Assumption 1 1 . 1 [ when N goes to infinity, one 
has 

T ATn T a.s. 
One may approximate t by t A Tjv if N is chosen such that Pp,{T > T/v) be small 



enough (the choice of N will be discussed in section 2.3 1 because the evolution of the 
process beyond Tj^ will have little impact on the law or the moments of the exit time. 
The rest of this section presents the two problems we are interested in and describes 
recursive methods to solve them. 



Definition 2.2. Let us introduce u* (x) for all x G U as the time for the flow starting 
from the point x to exit from U : 

u*{x) = inf {s > such that $(a;,s) ^ U} 

We now introduce some technical assumptions that will be in force throughout the 
whole paper. The first three ones will be crucial while the two last ones can be made 
without loss of generality. 

Assumption 2.3. The function u* is 

a. Lipschitz continuous, 

b. bounded by Cu* ■ 

Assumption 2.4. For all m e M , the set Um is convex. 

Assumption 2.5. For a > 0, let [/" — {x ^ E such that d{x,dU) < a}. There exists 
C > and /? > such that for all k e {0, N}, Pf,{Zk G f/") < Ca^. 

Remark 2.6. This technical condition can be checked in most of the applications. We 
will see, in the examples developed in Section [4[ how it can be derived quite generally 
when Zk has a bounded density. Moreover, it could be replaced by the following 
one, similar to an hypothesis introduced by M.H.A. Davis in (Section 24) and 
presented as quite general in applications : there exists e > such that for all x ^ U , 



10 



A. Brandejsky, B. de Saporta, F. Dufour 



QiW^x) = where U' ^ {x £ E such that d{x,dU) < e} i.e. for all k G {0, ...,N}, 

'Pf.iZk e U^) = 0. 

Assumption 2.7. The process cannot go back to U once it has left it i.e. Vz S 
C/",P.(Vt>0,Xtef/) = 0. 

Assumption 2.8. The function t* is bounded by Cf. 

In our discussion, Assumption |2 .7|does not imply any loss of generality and Assump- 



tion 2.8 stems from Assumption |2 .3| b . Indeed, if any of the two previous assumptions is 
not satisfied by the process (At)tgR+, we introduce the process killed at time r denoted 
{Xt)tm+ and defined by: 

{Xt for t <T, 
A for t > T. 

where A denotes a cemetery state. The state space of the killed process is E = UU{A} 



and Assumption 2.7 is fulfilled since the killed process remains in A after leaving U . 
In addition, i* , the deterministic exit time from E for the killed process equals u* that 
is bounded and Lipschitz continuous according to Assumption |2.3[ 

2.1. Distribution 

The first goal of this paper is to compute an approximation of the law of the exit 
time T. More precisely, we intend to approximate P^{t > s|t < Tjv) for s > 0. 

Our approach has a huge practical interest because we will see that, after some 
initial computations, any value of the survival function of t may be quickly obtained. 
More importantly, our approach is even flexible with respect to U in the sense that 
the survival function of the exit time r' from a new set U' d U will also be directly 



available (provided that Assumptions 2.3 to 2.7 are still fulfilled by [/') 



Definition 2.9. For all s > 0, define as follows the sequences {pk{s))k>Q, {qk)k>a Sind 
{rk{s))k>o 

{Pk{s) =Pp(r >s|r<Tfe), 
qk = Pp(T < Tfe), 
rfe(s) =Pp({r >s}n{Tfc <T<Tfc+i}). 
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Remark 2.10. The conditional probability pk{s) does not exist when — 0. We 
then choose to extend the sequence by setting pk{s) = 0. 

Our objective is to approximate Pn{s) where N represents the computation horizon. 
The following proposition provides a recursion for the sequence {pk)k<N, pointing out 
that Pat may be computed as soon as the sequences {qk)k<N and {rk)k<N-i are known. 



Proposition 2.11. Under Assumption 2.1 for all k s > 0, po(s) — and 
Pk+i{s) = 



Pk{s)qk+rk{s) 
Qk + l ' 

otherwise. 



if qk+i 7^ 



Proof. First, recall that Tq = so that one has pq = since the process starts in 
U according to Assumption |2.1[ Then, let A; G N such that qk+i 7^ and notice that 
{t < Tk+i} = {t< Tk} D{Tk<T< Tfc+i}, one has 



Pk+iis) 



P^{{r>s}n{T<Tk+i}) 
P,,(r < Tk+i) 

^ Pp({r > 5} n {r < Tk}) + P^({r > s} (1 {Tk < t < Tk+i}) 

Qk+l 

_ Pk{s)qk + rkjs) 

Qk+l 

showing the results. □ 

Now, before turning to computations, let us present the second problem we are 
interested in. 

2.2. Moments 

Our second goal is to approximate the moments of the exit time from U i.e. for 
all j e N, we are interested in Ep[t^ |t < Tjv]. This is a very classical problem 
and some results are already available. First, it is possible to use a Monte Carlo 
method and we will point out why the method we propose is more efficient and flexible. 
Furthermore, K. Helmes, S. Rohl and R.H. Stockbridge introduced in 9 a numerical 
method for computing the moments of the exit time based on linear programming. J.- 
B. Lasserre and T. Prieto-Rumeau improved this method in [10 by using semidefinite 
positivity moment conditions. These methods are quite efficient but they require an 
assumption related to the generator of the process which is generally not fulfilled by 
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the PDMP. The method we are introducing now is based on the use of the Markov 
chain {&n)neN = (^ri,Tri)neN associated to the continuous-time process {Xt)t£M+- 

Definition 2.12. For all j £ N, introduce the sequences {pk,j)k>Q aiid {rk.j)k>Q 
defined as follows 

Pk,j = [r^jr < Tfe] , 

rkj ^'Ef,[Tn{T,<r<n + ,}]- 

Our objective is to approximate pj^j where N still represents the computation horizon. 
Similarly to the previous section, the sequence {pk.j)k<N satisfies a recursion which 
parameters are the sequences {qk)k<N, previously introduced, and {rkj)k<N-i- 



Proposition 2.13. Under Assumption 2.1 one has for all k,j e N, poj = and 

otherwise. 



Pk+i,j = 



Proof. The proof is similar to the previous one (see proposition 2.11 1. □ 



Before turning to the approximation method itself, let us discuss the crucial question 
of the computation horizon. 

2.3. The computation horizon 

In this paragraph, let us study more precisely the construction of the process (Xt) 
in order to obtain some results concerning the jump times (Tk)keN. For this purpose, 
we introduce, in this section only two additional hypothesis. 

Assumption 2.14. The jump rate A is bounded by C\. 

Assumption 2.15. There exists e > such that for all x £ E, Q{x,Af) — 1 where 
— {x £ E such that t*(x) > e}. Roughly speaking, the jumps cannot send the 
process too close to the boundary of E. 

Assumption |2.14| is satisfied in a large majority of applications ; Assumption |2.15| is 
quite general too and was introduced by Davis in [S], section 24. 



Let (ri,.4, P) be a probability space on which is defined a sequence (nfc)fegN of 
independent random variables with uniform distribution on [0; 1]. Let x — (m,^) G E 
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and oj e f2 and let us focus on the construction of the trajectory t>Q} of the 

process starting from point x. Let 

L if t < 0, 

F(t,x)=<j exp(-/o A(m,*„(e,s))ds) ifO<t<r(a;), 
iit>t*{x). 



It is the survival function of the first jump time Ti. Define then its generalized inverse: 
*(u,a;) = 



ini{t > : F{t,x) < u}, 

+00 if the above set is empty. 



Let then S'i(w) = Ti(w) = *(ni(w),a;) and for all t < Ti(w), 

Xt{uj) = {m,<^m{^,t)). 

If Ti{lo) < +0O, choose Xt^ with distribution Q{.,^rn{i,Ti)). Assume the trajectory 
is constructed until time Tfc. If Tk{oj) < +oo, let 

5fe+i(a;) = *(nfe(a;),XTj, 
Tk+i{w)=Tk{oj) + Sk+i{co). 

If Tfe+i(cj) < +00, choose ^Tfc+i with distribution Q{.,^rnT^{S,T^, Sk+i))- The trajec- 
tory is finally constructed by induction. 

With the same notations as above, we may then state the following lemma : 

Lemma 2.16. Let H be a survival function such that for alltGM. and for all x G E, 
H{t) < F{t, x). There exists a sequence of independent random variables {Sk)keN with 
distribution H and such that 

yK e M, Viv G N, p^(Tjv < i^:) < P^(r,v < k). 

where Tjv = X^f=o -^fc- 

Proof. Let H be such a survival function and let "if be its generalized inverse i.e. 



inf{t > : H{t) < u}, 

+00 if the above set is empty. 
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The assumption made on H yields for all x E E, 'i'{u) < 'i>{u,x). Let for all fc £ N 
and for all w G 17, 

^feH = $(nfe(w)). 

Notice that we are using the same 11^ as in the definition of S^, allowing us to write 
that Sk < Sk a.s. and therefore Tk < Tk a.s. The result follows. □ 

Similarly to M.H.A. Davis (Section 33 in [^), we approximated r by r A since 
tKTm t as N +00 thanks to Assumption It is therefore necessary to choose N 
large enough such that Pfj_(Tj^ < r) be small. It is tough to estimate this probability for 
a general process because the links between r and the jump times are largely problem- 
dependant. For instance, the geometry of U can be very complex. Therefore, N will 
generally be estimated through simulations. Indeed, one may compute P^(T]^ < r) 
for some fixed N thanks to a Monte-Carlo method and increase the value of N until 
this probability becomes small enough. 

However, we introduce an other method to bound this probability that may prove 
useful in applications. First, notice that, for any K > 0, 

{Tjv < r} C {Tn <K}U{t> K}. 

This implies that 

P^.{Tn <t)< P^{Tn < K)+P^{t > K). 

This will prove especially useful whenever r is bounded, which happens quite often in 
application, because there exists then K such that P ^{t > K) — 0. On the contrary, 
when r is not bounded, it remains however sometimes possible to obtain K such that 
P^(t > K) be small. 

Example 1 : a crack propagation model We adapt here an example studied by J. 
Chiquet and N. Limnios in which models a crack propagation. Yt is a real-valued 
process representing the crack size satisfying : 

Yo>0 

Yt = AtYt for all t > 0. 

where At is a Markov process with state space where < a < jS . We are 

interested in the time r before the crack size reaches a critical size ye- Xt = {At,Yt) 
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is a PDMP, At representing the mode at time t. It is possible to bound the exit time 
by considering the slowest flow : one clearly has for all t > 0, Yt > Yoe"* o,^d, thus 

We now intend to bound Pp(T;v < K) for a fixed K > 0. Let 

1 if i < 0, 

H{t) = { e-^^* if < i < e, 
iit>e. 

Distribution H represents, roughly speaking, the worst distribution of the inter- 
jump times Sk in the sense that it is the one that imphes the most frequent jumps. 
Indeed, denote Fk the survival function of Sk, one has H < Fk for all k E N. Therefore, 
Lemma 2.16 provides a random variable Tjv — '^k=o where Sk are independent and 
have survival function H such that 

P^{Tn <K)< P^(TW < K). 

We now bound P^(r;v < K). Standard computations yield E^[T;v] = Nm and 
V^JTat] = Na^ where : 



Assume now that N is such that Nm > K and notice that 

Tchebychev inequality yields : 

F,{T^<K)<- 

[Nm — K) 

and that the right-hand side term goes to zero when N goes to infinity. 



Finally, when r is bounded with a high probability and when Assumptions |2.14| and 
2.15 are fulfilled, we are able to choose N a priori such that V ^{Tjsi < r) be small. 
These conditions are satisfied in a large class of applications. 
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3. Approximation scheme 

3.1. The quantization algorithm 

First of all, let us describe the quantization procedure for a random variable and 

recall some important properties that will be used in the sequel. There exists an 

extensive literature on quantization methods for random variables and processes. We 

do not pretend to present here an exhaustive panorama of these methods. However, the 

interested reader may for instance, consult the following works [T1[7J[TT] and references 

therein. Consider X an M'- valued random variable such that \\X\\ < oo where \\X\\ 

II lip II lip 

denotes the Lp-nom of X: \\x\\^ = (^E[|X|p]) . 

Let Khe a fixed integer, the optimal Lp-quantization of the random variable X consists 
in finding the best possible Lp-approximation of X by a random vector X taking at 
most K values: X S {x^, . . . ,x^}. This procedure consists in the following two steps: 

1. Find a finite weig hted grid F C M'' with F = {a;\ . . . , x^}. 

2. Set X = X^ where X^ — projr{X) with pr denotes the closest neighbour 
projection on F. 

The asymptotic properties of the Lp-quantization are given by the following result, see 

e.g. m- 

Theorem 3.1. Jf E[|X|?'+'''] < +oo for some 77 > then one has 

lim min \\X ~ X^\\p = J, 

K^oo \r\<K P ^ 

where the law of X is Px{du) = h{u)Xq{du) + v with v L Xd, Jp^d a constant and Xq 
the Lebesgue measure in M'^. 

Remark that X needs to have finite moments up to the order p+rj to ensure the above 
convergence. There exists a similar procedure for the optimal quantization of a Markov 
chain {X/jjfegN- There are two approaches to provide the quantized approximation of 
a Markov chain. The first one, based on the quantization at each time k of the random 
variable Xk is called the marginal quantization. The second one that enhances the 
preservation of the Markov property is called Markovian quantization. Remark that for 
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the latter, the quantized Markov process is not homogeneous. These two methods are 
described in details in [111 section 3]. In this work, we used the marginal quantization 
approach for simplicity reasons. 

Our approximation methods are based on the quantization of the underlying discrete 
time Markov chain {Qk)k<N = iZk,Tk)k<N- The quantization algorithm provides for 
each time step < fc < a finite grid Ffc of _E x IR+ as well as the transition matrices 
{Qk)a<k<N-i from Tk to r^^+i. Let p > 1 such that for all k < N, Zk and Tk have 
finite moments at least up to order p and let projri. be the nearest-neighbor projection 
from E X M+ onto Tk- The quantized process {Qk)k<N = {Zk,Tk)k<N with value for 
each k in the finite grid Vk oi E x is then defined by 

{Zk,Tk) ^ projr^{Zk,Tk). (3) 

In practice, we begin with the computation of the quantization grids which merely 
requires to be able to simulate the process. These grids are only computed once and for 
all and may be stored off-line. Our schemes are then based on the following simple idea: 
we replace the process by its quantized approximation within the different recursions. 
The results are obtained in a very simple way since the quantized process has finite 
state space. 

Remark 3.2. In addition, we recall a technical property of the quantization algorithm 
proved by C. Bouton and G. Pages in ^'3' ; the quantized process evolves within the 
convex hull of the support of the law of the original process. Therefore, and it will be 



required below. Assumption 2.4 yields that, if Zk a.s. for some fc g {0, N} then 
Zk & U a.s. 

3.2. Approximation scheme of the distribution and proof of convergence 



We already noticed in Proposition 2.11 that Pn{s) = P^iT > s\t < Tjy) may be 



computed as soon as the sequences {qk)k<N and {rk)k<N-i are known. Therefore, we 
will find expressions of these sequences depending on the Markov chain {Zk,Tk)k<N 
that we will replace by the quantized process {Zk,Tk)k<N in order to define their 
quantized approximations {qk)k<N and {rk)k<N-i- 

First, notice that {Tk < t} ^ {Zk G U} and {r < Tk} = {Zk i U} thanks to 
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Assumption [2?7| Moreover, on {Z^ G U, Z^+i <^U}, one has r = (T^, + u*{Zk)) A T^+i 
a.s. where u*{x) is the deterministic exit time from U starting from the point x (see 



Definition 2.2|), and one has : 

rfe(s) = 'E^j\l{^Tk+u'{Zk))fyTk + i>s}'^u{Zk)lu''{Zk+l)\- 

The above equations are crucial in our discussion and, from now on, we will use them 
without referring to Assumption |2.7[ 

Before turning to the approximation scheme itself, let us state some properties of the 
sequence {qk)k<N that will be important in the following proofs. Indeed, the sequence 
{qk)k increases since {t < T^} C {t < Tk+i} for all fc < — 1. Moreover, note that 



qo — and lim„_).-|_oo 9ji = 1 thanks to Assumption 2.1 Therefore, there exists an 
index denoted fc > 1 such that 

• for all k < k, one has — 0, 

• for all k > k, one has qk > 0. 

We denote q = qj. the first positive value of the sequence so that qk > q for all fc > A;. 
One obtains then the following definition. 

Definition 3.3. Let 

fc = inf {fc > such that qk > 0} , 

i.e. q is the first strictly positive value of the sequence {qk)ke{o,...,N}- 

We now naturally define the quantized approximations of the previous sequences. 

Definition 3.4. For all s > 0, define the sequences {qk)ke{o,...,N} and {rk)k<£{o,...,N-i} 
by: 

' qk =E^[l[/c(Zfe)], 
It is important to notice that both qk and ffc(s) may be computed easily from the 
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quantization algorithm. Indeed, one has : 

qk= = 

e = {z,t)(iVk 



zeU z' 4U 



Recall from Proposition 2.11 that the sequence (j)k)k<N satisfies a recursion depend- 
ing on two parameters: {qk)k<N a-nd {rk)k<N-ii that we are now able to approximate. 
Hence, replacing them by their quantized approximations within the same recursion 
leads to a new sequence denoted {pk)k<N- The rest of this section is dedicated to the 
proof of the convergence of {pk)k<N towards {pk)k<N- This convergence is far from 
being trivial because on the one hand, the definitions of the sequences {qk)k<N and 
i'>'k)k<N-i contain many indicator functions that are not Lipschitz continuous and on 
the other hand, the recursive function giving Pk+i from pk, qk, qk+i and is not 
Lipschitz continuous either. 

Definition 3.5. For all s > and for all k e {0, N - 1}, let po{s) = and 



Pk+iis) 




ifqk+l^O 
otherwise. 



(5) 



The two following propositions will be necessary to prove the convergence of the ap- 
proximation scheme. They respectively state the convergence of {qk)k<N and {rk)k<N-i 
towards {qk)k<N and {rk)k<N-^i- 



Proposition 3.6. Under Assumptions 2.5 and 2.7 for all k e {0, N}, qk converges 
towards q^ when the quantization error — QkWp goes to zero. More precisely, the 
error is bounded by 



\qk-qk\ < Cp+f 




\Zk-Z, 



I p+3 



k\\p 



where C and (3 are defined in Assumption 2.5 
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Proof. For all k e {0, iV}, equation Q yields 

\qk-qk\ = \Ef,[lu{Zk)-lu{Zk)]\- 

The difference of the indicator functions is non zero if and only if and Zk are 
on either side of dll. Therefore, in this case, for all a > 0, if \Zk — Zk\ < a, then 
d{Zk,dU) < a. Hence, either — Zk\ > a or Z^ <E U". Markov inequality and 
Assumption |2.5| yield : 

^^\lu{Zk) - luiZk)\ <P^{\Zk - Zk\ > a) +P^{Zk e U") 



< 



Zk\ 



aP 



This bound reaches a minimum when a = ( ''^^ | and the result follows. □ 



Proposition 3.7. Under Assumptions 2.3 a, 2.5 and 2.7 for all k G {0, N — 1} and 

for almost every s > w.r.t. the Lebesgue measure on M, 

when the quantization errors \\Qi — Qi\\p for I e {fc, fc + 1} goes to zero. 
Proof. Let k G {0, ...,N - 1} and s > 0, equation Q yields 

\rk{s)-rk{s)\ <A + B. 

where 



A = 
B = 



E„ 



^{{Tk+u*{Zk))ATk+i>s} 



^{(.¥.+uHZ,))A¥.^r>s})'^u{Zknu^{Zk+l) 



^^^^^(luiZk)lu4Zk+i) - lu{Zk)lu4Zk+l) 



{(Tfc+«'(Zfc))ATfc^ 

In the A term, we crudely bound lu{Zk) and lijc(^Zk+i) by 1 and turn to the difference 
of the two indicator functions. This difference is non zero if and only if (T^ + u*{Zk)) A 
Tk+i and (Tk + u*{Zk)) A Jfc+i are on either side of s yielding that they both belong 
to [s — 77; s + rj] where r] — \(Tk + u*{Zk)) A Tfc+i — {Tk + u*{Zk)) A 7fc+i I , one has then: 



|l{(Tfc+«-(Zfc))ATfc+i>s} - l{(Tfc+«-(Zfc))Affc+i>s}l - ^{\iTk+u-{Zk))^n + i-s\<T,} 



so that 



A<P^{\{Tk + u*{Zk))ATk+i-s\ <f]) 
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The following discussion consists in noticing that either i] is small and so is the 
probability that (T^ + A Tfe+i belongs to the interval [s — ij, s + r]], or 77 is large 

but this happens with a small probability too when the quantization error goes to zero. 
For all a > 0, one has 



A < {\(Tk + u*{Zk)) A Tfc+i - s| < r/, ?7 < a) + (77 > a) 
< (KTfc + u*{Zk)) A Tfe+i - s| < a) + P^ (77 > a) 



< \ipk{s + a) - (fk{s - a)\ + 



where (pk denotes the distribution function of [T^ + u*{Zk)) A Tk+i- Let e > and 
assume that s is not an atom of this distribution so that there exists ai > such 



that \ipkis + ai) — ipk{s — ai)\ < e. Besides, thanks to Assumption 2.3 a stating the 



Lipschitz continuity of u* , one has rj < \Tk — Tk\ + [u*]\Zk — Zk\ + \Tk+i — Tk+\\ 
Moreover, since the quantization error goes to 0, one may assume that ||?7||p < aie? 
Setting a = ai in the previous computations yields 



A < \ipk{s + ai) - ipkis - Q!i)| + 



1^/11 



< 2e. 



Notice that the set of the atoms of the distribution function of (T^ + u*{Zk)) A Tk+i is 
at most countable so that the previous discussion is true for almost every s > w.r.t. 
the Lebesgue measure. Let us now bound the B term: 



B < EJlu{Zk)lu4Zk+i) - luiZk)lu4Zk+i)\ 



<E^ \lu^-{Zk+i)\lu{Zk)-lu{Zk) 
< \qk- qk \ + \qk+i - qk+i\ 

that goes to zero thanks to Proposition |3.6| 



^u{Zk)\^U''{Zk+i) — ^u^iZk+i)] 



□ 



The convergence of the approximation scheme of the distribution of the exit time is 
now a straightforward consequence of the following proposition. 



2.1 



2.4 



2.5 



and 2.7 hold. Let {(7k) 



k<N-l 



Proposition 3.8. We assume Assumptions 
and {ak)k<N-i be two sequences of [0, l]-vaJued real numbers. Let {TTk)o<k<N a,nd 
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iTfk)o<k<N defined as follows, ttq = ttq = and 




if Qh+ij^O 
otherwise. 



ifqk+i^O 
otherwise. 

For < k < N , if the quantization error is such that for all I < k 
then 




2 / 

a V 



2(7r''«P + 1) 



Qk - Qk 



where 7r''"P = maxo<fe<7v T^k- 



Proof. The difRculty of this proof hes in the fact that the recursive function giving 
TTfe+i from TTfc, qk, qk+i and ak is not Lipschitz continuous because of the division by 
qk+i- To overcome this drawback, we will use the strictly positive lower bound for q^ 
described earlier. Indeed, recall from Definition 3.3 that there exists a step k such that 



Qfc > 9 > for all k > k and qk = for all fc < fc. What is more, a similar bound will 
be derived for the quantized values qk thanks to the convergence of qk towards qk . 

We now prove by induction that tt^ converges towards tt^.. First, one has ttq = ttq = 
0. Then, let ke {!,..., TV}. 



If A; < fc, then qk = and Assumption 2.4 yields that qk = too. Indeed, qk = 
means that Zk G U a.s. Since J7 is a convex set. Remark 3.2 implies that Zk G U a.s. 
too. In other words, qk ~ 0. Finally, from the definitions, one has tt^ = tt^ = 0. 

If fc > fc, then qk > q > 0. In order to bound the error between tt^ and tt^, it is 
indeed necessary to have a strictly positive lower bound for qk because of the division 
by qk within the recursion. Now we need to obtain the same kind of bound for qk. 



This can be achieved thanks to Proposition 3.6 giving the convergence of qk towards 
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q^. Indeed, assume from now on that the number of points in the quantization grids 
is large enough such that the quantization error is sufficiently small to ensure that for 
all j — k,...,N, \qj — Qjl < ^q. Hence, the required lower bound is qk > -^q > 0. 
Therefore, 



qk qk 



^T^k-l I - I , I - I , I - I 

<^^|9fe-i - qk-i \ + -^^\nk^i - TTfc-il + ^|tTfc_i - ak-i\ 

qk qk qk 



\qk - qk\ 



<- 



qk 



\qk-i - qk-i\ 



qkqk 

1 , ^ , 1 , ^ , 
qk Qk 



+ (tt^-p + 1) 



\qk - qk\ 

qkqk 



<- (TT^^'l^fc-l - qk-l\ + \TTk-l - ^ffc-ll 

q 



1) 



\qk ~ qk\ 



where tt* 



maXo<fc<Ar TTfc. 



□ 



Remark 3.9. Notice that a bound for the rate of convergence ofifk towards TVk may 
he obtained as soon as a bound for the rate of convergence of ak towards Uk and an 
upper bound for the sequence {TTk)o<k<N are available. 

Eventually, let us state one of our main results, namely the convergence of the 
approximation scheme of the distribution of the exit time: 



Theorem 3.10. Under Assumptions \2J\ \Z3\ a, \2J\ and \27^ for all k e {0, ...,N} 
and for almost every s > w.r. t. the Lebesgue measure on M, 



Pk{s) Pk{s) 



when the quantization errors ||0j — Oj||p for j G {0, fc} go to zero. 

Proof. Let s > such that (ffc(s))fc converges towards {rk{s))k and apply Propo- 



sition 



3.8 with (crfc)fe = (rfc(s))fe and {ak)k = (rfc(s))fe so that (TTfc)* 



(7ffe)fc = (Pfe(s))fc. Finally, notice that (pfc(s))fc is bounded by 1. 



{Pk{s))k and 
□ 
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Remark 3.11. It may be useful to notice tliat, although it will be crucial in the 



moments approximation scheme, the boundedness condition on u* (Assumption 2.3 b) 
was unnecessary in this section. Hence, the distribution approximation can be achieved 
without this hypothesis. 

Eventually, we obtain an easily computable approximation of the survival function 
of the exit time. Let us now consider its moments. Of course they may be derived 
from the distribution but we present in the following section a method to approximate 
them directly. An important advantage of this method will be to provide a bound for 
the rate of convergence. 

3.3. Approximation scheme of the moments and rate of convergence 

Similarly to the distribution, the moments may be approximated thanks to the 
quantization of the process {Qk)k<N — {Zk,Tk)k<N ■ However, it is important to stress 
the fact that we will be able to derive a rate of convergence for our approximation 
scheme. One may notice from Proposition |2.13| that, similarly to the case of the 
distribution, pisi,j — Ep[T^ |r < T/v] may be computed as soon as the sequences {qk)k<N 
and {rkj)k<N-i ai'e known. The first one has already been approximated in the 
previous section but we still need to find an expression of the second one depending on 
the Markov chain (Zk, Tk)k to define its quantized approximation {jk,j)k<N-i- Thanks 
to Assumption |2.7| the same arguments give 



{{Tk + u*{Zk)) hTk+iyiu{Zk)lu^{Zk+i) 



(6) 



So that we may now naturally define the quantized approximation of the sequences 

{rk,j)k<N-i and {pk,j)k<N- 

Definition 3.12. For all j G N, define the sequence {rk.j)ke{o....,N-i}' 



TkJ = Ep 



((ffe + u*iZk)) A fk+iY lu{Zk)lu4Zk+i) 



and the sequence {pk,j)k€{o,...,N} by pqj — and 




if Qk+i ^ 
otherwise. 
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As for qk and ffe(s) defined in the previous section, r^ j may be computed easily from 
the quantization algorithm. Indeed, one has : 



E 

zee/ 



E 

{z',t') e Tfc+i 



(t + «*(z))At' v{Qk = e)Qk{e;e'). 



The following proposition proves the convergence of ffc towards r^^ 



Proposition 3.13. fJiider Assumptions, \2^a, ^^\Z7\and^^ for all k e {0, N ~ 

1} and for all j G N, j converges towards r^.j when the quantization errors \\Qi — Qi\\p 
for I e {k, k + 1} go to zero. More precisely, the error is bounded by 

\rk.j - rk,j\ < j{{k + \)Ct'y'\\\Tk - ffellp + [u*]\\Zu - ZuWp + \\Tk+i - n+i\\p) 
+ ((fc + l)Ct'Y{\qk -qk\ + kfe+l - ?fe+l|). 

Proof. Let k £ {0, N - 1} and j e N, one has : 

\rk,j-rk,j\ <A + B. 

where 
A= E„ 



((Tfe + u*{Zk)) A Tfe+i)' - ((Tfe + u*{Zk)) A Tfc+i)' j ly(Zfe)l[/c(Zfc+i) 
in + u*{Zk)) An+i)' (lu{Zk)lu4Zk+i) - lu{Zk)lu4Zk+i) 



E„ 



Assumption |2.8| yields that the inter jump times Si are a.s. bounded by Ct* so that 



T, < iCf a.s. and (T, + u*{Zi)) A T^+i < (i + l)Ct. a.s.. By using Remark pO) these 
bounds are equally true for the quantized process Ti < iCt* and (Ti + u* {Zi)) A T^+i < 
T.+i <{i + \)Ct. a.s. 

Let us first consider the term A, we crudely bound the indicator functions by 1. 
Moreover, denote 77 = \{Tk + u*{Zk)) A Tk+i - {Tk + u*{Zk)) A Tk+i\ and notice 
that the function x — >■ a;-' is Lipschitz continuous on any set [0, M] with Lipschitz 
constant jIv'P~^. 



A<E„ 



]{{k+i)Ct.y \ 



< 
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and thanks to Assumption |2.3[ a stating the Lipschitz continuity of u* , one has 
A < j{{k + \)Ct^Y'U\\Tk -Tu\\p+ [u*]\\Zk - Zfcllp + - ffc+illp 



Moreover, the term B is bounded by: 

< {{k + l)Ct.)-'(|gfe - gfcl + \qk+i - qk+i\). 
We conclude thanks to Proposition |3.6[ 



□ 



We may now state the other important results of our paper namely the convergence 
of the approximation scheme of the moments of the exit time with a bound for the rate 
of convergence. 



Theorem 3.14. Under Assumptions [2^a, [2^ |2j and for all k e 
{0, ■■■,N} and for all j G N, pk,j converges towards pk.j when the quantization errors 
\\Qj — QjWp for j e {0, k} go to zero. 

More precisely, if the quantization error is such that for all I < k 



Cp+1 



\Zi~Zi\\^' 



1 



then 



\Pk,j -Pk,j\ < -~ {{NCt'Y\qk-i - qk-i \ + bfc-ij -Pk-i.j \ + \rk-i,j -^^fe-ijl) 



2{{NCt>y + l) 



qk - qk\- 



Remark 3.15. The rate of convergence depends on the quantity q whose exact value 
might be unknown in some complex applications. In that case, it may still be approx- 
imated through Monte-Carlo simulations (see examples in Section [4]). Nevertheless, 
Theorems \3.10\ and \3.14\ prove the convergence of our approximation schemes regardless 
of the value of q. 



Proof. Let j G N and apply Proposition 3.8 with {<7k)k = {fk.j)k and {^k)k = {^k,j)k 



such that (7rfc)fc = {pk,j)k and {T:k)k = {Pk,j)k- Finally, according to Remark 3.9 a 
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bound for the rate of convergence is obtained since the sequence {pk,j)o<k<N is bounded 
by: 



Hence, the result. 



Tl\r<Tk 



< [{kCt>y\T < n] < {kCt,y < {NCt>y 



□ 



4. Examples and numerical results 

4.1. A Poisson process 

Let Nt be a Poisson process with parameter A = 1 and let Yt — t + Nt- (lt)t>o is 
a PDMP with state space E — R ; inter jump times Sk have independent exponential 
distribution with parameter A = 1 ; the flow is defined on (M+)^ by ^{x,t) — x + t ; 
and finally, the post-jump locations satisfy : Va; G -B, Q{{x + 1}, x) — 1. An example 
of trajectory of the process is represented in figure [T] We are interested in the exit 
time problem for the process {Yt)t>o. The study of this process is especially interesting 
because it is possible to compute the exact value of its distribution function in order 
to compare it with the numerical value given by our approximation scheme. 

20- 
18- 
16- 
14- 
12- 
10- 

8- 

6- 

4 - 

2- 

0- 


Figure 1: A trajectory of the process (Yt) drawn until the 10*'' jump time. 
Let us turn now to the numerical simulations. Let 6 = 10 i.e. U ~] — oo, 10[. We 




1 2 3 4 5 6 7 8 5 10 
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may choose iV = 10 since — T, 



N 



T-i 



N 



N > N. Besides, it is clear that 



for all y g] — oo, 10[, u*{y) = 10 — y. Assumptions 2.3 and 2.4 are clearly satisfied and 
so is Assumption |2.5| thanks to the following lemma. 



Lemma 4.1. For all a > and for all k e {0, N}, 

(Zk e C/") < 2a. 

Proof. Since Zq = a.s., (Zq G C/") = (Zq G [10 - a, 10 + a]) = l{„>io} < 
j^a < 2a. 

Let now k G {1, N}. Denote f'y{k,i) the density of the distribution 7(fc, 1) and denote 
^ (fc-i)! (^7^)*^ ^ its bound. Since has distribution 7(fc, 1), Zk = k + Tk has 
density fz^i') — f-y{k,i)(' ^ ^) that is also bounded by Cfc. Eventually, one has: 

P,. {Zk G C/") = P^ (Zk G [10 - a, 10 + a]) < 2Cka < 2a. 

Indeed, the sequence {Ck)k decreases so that for all k G {1, ..-^N}, Ck < Ci — 1. □ 



Moreover, Assumption |2.7| is satisfied since the process increases but Assumption 
2.8 is not, because t*{x) = +oo for all x <E E. However, as pointed out in Section [2| 
this can be solved by considering the process killed at time r. 

The mean exit time Table [l] presents the simulations results for the approximation of 
the mean exit time. It includes for different number of points in the quantization grids 
the value of pn,i which approximates the mean exit time. A reference value is obtained 



thanks to Monte Carlo method (10^ simulations): E[tio] 



Monte Carlo 



5.125. 



Points in the quantization grids 




relative error to 5.125 


20 points 


5.050 


1.46 % 


50 points 


5.096 


0.56 % 


100 points 


5.095 


0.58 % 


200 points 


5.118 


0.13 % 


300 points 


5.128 


0.06 % 


500 points 


5.123 


0.03 % 



Table 1: Simulations results for the mean exit time 
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The second moment We present the results of the approximation of the second moment 
in Table |2] Our Monte Carol refer ence value (10*' simulations) is E[T2Q]Monto Carlo — 
27.5. 



Points in the quantization grids 


PN,2 


relative error to 27.5 


20 points 


26.66 


3.05 % 


50 points 


27.20 


1.11 % 


100 points 


27.21 


1.05 % 


200 points 


27.43 


0.25 % 


300 points 


27.54 


0.13 % 


500 points 


27.49 


0.03 % 



Table 2: Simulations results for the second moment 



For the first and second moment, the empirical convergence rate is presented on 
Figure [2] It is estimated through a regression model as —1.23 for the first moment and 
— 1.39 for the second moment. Remark that there are roughly the same order as the 



rate of convergence of the optimal quantizer (see Theorem 3.1 1 as here the dimension 
is 1. 




Figure 2: Logarithm of the error w.r.t. the logarithm of the number of points in the 
quantization grids for the first and second moment of the Poisson process. 
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The exit time distribution As announced earlier, one can obtain the exact value of the 
survival function of the exit time. 

Proposition 4.2. Denote fl{.) the floor function. For all s,b E M+, one has : 

P[Tfi(^b-s)+i > s] foralls<b, 



P[n >s] = 







otherwise. 



Remark 4.3. Notice that has distribution ^{k, 1) so that the right-hand side term 
in the above proposition can be computed easily. 

Proof. Let s > 0. Notice that Yg > s, thus < s a.s. when s > b. Assume now 
that s < b, one has : 

P[n > s] = P[n <b] = P[N, <b-s] = P[N, < fl{b - s)] = P[r^,(,_,)+i > s]. 

Hence, the result. □ 

Figure[3]represents both the exact survival function of the exit time and its quantized 
approximation. Table [3] contains the empirical error between the two functions. For 
the survival function, the empirical convergence rate is presented on Figure [4] It is 
estimated through a regression model as —1.05. Remark that it is roughly the same 



order as the rate of convergence of the optimal quantizer (see Theorem 3.11 as here the 
dimension is 1. 



Points in the quantization grids 


maxs \pn{s) -Pn{s)\ 


20 points 


0.090 


50 points 


0.077 


100 points 


0.057 


200 points 


0.011 


300 points 


0.007 


500 points 


0.005 



Table 3: Simulations results for the distribution 



Remark 4.4. We already insisted on the fact that our approach is flexible w.r.t. U . 
In this example, one could obtain very quickly the mean exit time or the exit time 
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Figure 3: Survival function of no and its quantized approximation with 500 points in the 
quantization grids. The functions appear indistinguishable. 



10 




1 : 
10 10 



Figure 4: Logarithm of the error w.r.t. the logarithm of the number of points in the 
quantization grids for the survival function of the Poisson process. 
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distribution for a different set U' ~] — oo, b'] for any < &' < 5 = 10. Indeed, 
P(rb/ > Tio) = so that it is not necessary to compute new quantization grids. 

Remark 4.5. Recall that the value ofTj, may be obtained from since = ~ k 
so that it is sufRcient to quantize the process {Zk)k<N instead of {Zk,Tk)k<N- The 
reduction of the dimension of the process that has to be quantized results in an 
improvement of the convergence rate and it appears that the approximations presented 
in the previous tables converge indeed very quickly. 

Convergence rate for the exit time distribution One may notice from the proof of 
Proposition |3.7| that a bound for the rate of convergence for the exit time distribution 
may be obtained as soon as for all k E {0, — 1}, the survival function of {Tk + 
u*{Zk)) A Tk+i denoted tpk is piecewise Lipschitz continuous. Although it is tough to 
state general assumptions under which this is true, the following proposition proves 
that the condition is fulfilled in our example. 

Proposition 4.6. Forallk e {0, A^— 1}, the survival function ipk of {Tk + u* (Zk)) A 
Tk+i is Lipschitz continuous on ] — oo; b^k[ and on ]6— fc; +oo[ with Lipschitz constant 

K] < 1. 

Proof. Let k = and s > 0, one has: 

Ms) = P^((ro + u*(Zo))Ari >s) 

= Pf,{b/\Ti>s) 
= l{fc>,}P^(Ti > s) 

= l{fc>s}e^'' since Ti has exponential distribution with parameter 1. 

Therefore, the function (po is worth zero on [6;+oo[ and is Lipschitz continuous with 
Lipschitz constant 1 on ]0;b[. 

Let fc > 1, s > and remember that the random variables (S'j)j>o are independent 
and all have exponential distribution with parameter 1 so that, in particular, Tk and 
5^-1-1 are independent and Tk has distribution 7(fc, 1). Moreover, recall that Zk — k+Tk 
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and that u*(x) = b — x. 

(^fe(s) = P^((Tfe + u*{Zk)) A Tfe+i > s) 

^{{t+(b~k~t))Au>s} f-,(k.i){t) f-t(k+i,i){u)dtdu 

where denotes the density function of the distribution 7(j, 1) for j G {k, k + 1}. 

Let s' > s > 0, one has: 



\(pk{s') - ipk{s)\ < / |l{(&-fc)An>s'} - l{(b-fc)AM>s}| /7(fe.l)(^)/7(/£ + l,l)(")'^^'^'^ 

J(R+)2 

< / l{(6-fe)AnG]s;s']}/7(/£,l)(^)/7(fc+l,l)(")o?*'^U 

J(R+)2 

< / (l{6-/£e]s;s']} + l{ue]s:s']}) /7(fe,l)(^)/7(*:+l.l)(^)'^*'^'" 

J(R+)2 

1 /fc^ 



< l{b-fce[s;s']} + |s - s| Since = — j < 1. 

Eventually if s and s' both belong to ]0; b — k[ or if they both belong to ]b — k; +oo[, 
one has \(pk{s') — ipk{s)\ < \s' — s\. The result follows. □ 



Consequently, in this example, we are now able to state a bound for the rate 
of convergence of the exit time distribution approximation scheme. The following 
proposition is therefore an improvement over Proposition |3.7| and Theorem |3.10[ 

Proposition 4.7. For all k E {0, iV— 1}, iet s > and assume that the quantization 
error be small enough to ensure that 

{\\Tk-fk\\p + \\Zk-Zk\\p+\\n+i-fk+i\\p)^ <\b-k-sl 

one has then 

kfe(s)-ffc(s)| <2(^)^ fl + l) (\\n-fk\\p + \\Zk-Zk\\p+\\Tk+i-fk+i\\p 



.27 \p 

+ \qk -qk\ + \qk+i - qk+i\- 

Moreover, for all k G {0, ...,N}, if the quantization error is such that for all I < k 
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one has then 

,2 

q 



\Pk{s) -Pk{s)\ <|(|gfc-i - qk-i\ + bfc-i(s) -Pk-i{s)\ + |rfc_i(s) - ffc_i(s)|^ 
4 , 

+ ^Wk - qk\- 

q 

Proof. The proof derives directly from the proofs of Proposition |3.7| and Theorem 



3.10 Simply notice that the A term may be bounded thanks to the piecewise Lipschitz 
continuity of the functions (pk on ] — oo; b — k[ and on ]b — k; +oo[. Let s > 0, s ^ b — k, 
and let a > such that b ~ k ^ [s ~ a; s + a] i.e. a < |6 — fc — s|, one has 

^ < \ fk{s + a) — (pk{s — a) H from the proof of Proposition 13.71 

< 2[(pk]a+ —f- 

■ Notice besides that Ypk\ = 1 and 
\u*\ = 1. □ 

Remark 4.8. We can calculate the exact value of q that is the first nonnegative value 
of the sequence (P^(Zfe ^ C/))^,. One has q = P^(Zi ^] - oo; 10[) = P,,(Ti > 9) = e^^ 
because Ti has exponential distribution with parameter 1. 

4.2. A corrosion model 

Let us consider a structure of aluminium corroded successively into 3 different 
environments. Corrosion is prevented by some protection until a random time 7 when 
corrosion starts. Then, in each environment i G {1; 2; 3}, the loss of thickness satisfies: 



lit) ^ p, (t - "f + 1], (e - 1)) l{f>^} 



where pi is the corrosion rate (pi has a uniform distribution on an interval that depends 
on the environment i) and ?/; is a constant transition time. The structure goes from 
environment 1 to environment 2, then from 2 to 3, from 3 to 1 and so on. It remains 
in environment i for a time Ti which has exponential distribution with parameter A.^. 
When the loss of thickness reaches 0.2 mm, the piece is said to be unusable, this will 
be the exit criterion. Table |4] gives the values of the different parameters. 



The loss of thickness will be represented by a PDMP whose modes are the different 
environments. Let then M — {{i,j) : i € {l,2,3},j E {0,1}}. For m — {i,j) G M, 
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environment 1 


environment 2 


environment 3 




(17520)-! 


(131400)-! 


(8760)-! 


m (h) 


30000 


200000 


40000 




[10-6,10-5] 


[10-^10-6] 


[10-6,10-5] 


1 (h) 


WeibuU distribution with a = 2.5 and /3 = 11800 



Table 4: Numerical values of the parameters of the corrosion model 

i represents the environment and j is worth 1 if the protection 7 is still active and 
otherwise. For each m e Af, let = and for ^ e E„i, ^ represents the family 
(d, s, p, 7) where d is the corroded thickness and s is the time since the last jump. The 
set Um will therefore be for all m £ M, U„i =] — 00; 0.2] x M?' . This set is convex so 

( 



that Assumption 2.4 is satisfied. Finally, the flow in mode m — (i, j) is : 



*(i,0)( 



d 

S 

p 

/o\ 

s 

p 

V ^ J 



d+ dyn{t+ S) - dm{s) 
t+ S 



\ 



p 







t + s 
P 

\ (7-^)l{7>t} j 

The parameters d and 7 evolve continuously between the jumps but p is chosen 
independently after each jump and is constant along the flow. 



Let us consider the approximation of the distribution and of the mean exit time. 
Concerning the first moment, one may notice that E^[t] = £^[7] + Ep[T'] where 7 has 
WeibuU distribution and r' represents the exit time in the case of a process without 
initial protection against corrosion (i.e. 7 = 0). Therefore, it is sufficient to check 
whether r' satisfles the required assumptions. Hence, let 7 = and notice that u* 
is then bounded since p > 10^^ and ry < 200000 so that d„,it) > IQ-'^^t - 200000) 
and eventually u* < 0.2 x 10^ + 200000 = 2.2 x 10^ h. Denote by C„. this bound. 



Concerning the distribution. Assumption 2.3 b (the boundedness condition on u*) is 
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not required according to Remark |3.11| Moreover, it is easy to notice from the proofs 



of Propositions 3.7 and 3.13 that Assumption |2.3|a (the Lipschitz continuity condition 



on u*) becomes useless in this example thanks to Lemma 4.10 Assumption 2.5 follows 



from Lemma [4. 9| below. Eventually, Assumption |2.7| is satisfied but Assumption |2.8| is 
not. However, considering the process killed at time t solves this issue. 

Lemma 4.9. For all a > and for all k e {0, ■■■,N}, 

Ft^iZk € < 5a. 

Proof. For notational convenience, introduce Mk, Dk, Rk and Gk the values of m, 
d, p and 7 after the /c-th jump so that — (Mk, Dk, Rk, Gk)- Notice now that 

Pt.{Zk e [/") = P^{\Dk - 0.2| < a). 

We therefore study more precisely the law of Dk- Let K = inf{fc > such that Gk = 
0}, K is the jump happening at the end of the protection against the corrosion. 
Eventually, denote F{s) = s + rj{er^ — l). One has then 

Dk = Q for fc < K, 

Dk = Dk^i + RkF{Sk) for k>K. 

Let us now prove that for all fc, the random variable RkF{Sk) has a bounded density. 
Recall that Rk has a uniform distribution on [ak;hk] C [10"^; 10~^] and Sk has an 
exponential distribution with parameter A^. Let now ft, be a real bounded measurable 
function. 



E^[/i(i?feF(5fc))] 



Introduce the following transformation 



h{pF{s))-^ 

afc Ok - ak 



Xke-^"' 



dpds 



u = p 

V — pF(s) 



whose Jacobian is worth ^{F ^)'(^) so that 



B^[hiRkF{Sk))] 



+ 00 



h{v) 



Xke-^-^-'(iHF-'y{l) 

{bk - ak)u 



du dv. 
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Hence, we obtain the density of the random variable RkF{Sk) and integration by parts 
yields 



Jau {h - ak)u bk - ak Ja^ 



u X ^ —du 



ue 



-AfcF-i(^) 



bk - ak 

Finally, the density of the random variable RkF{Sk) is bounded by 



du 



< 



ak 



5fc - Ofe 



1 < 



2bk 



< 2. 



5fc - flfe 



Let j £ N, wc now study the distribution of the random variables {Dk)keN conditionally 
to the event {K = j}. An induction argument provides that, conditionally to the event 
{K = j}, the random variable Dk has distribution Sg for k < j and has a density ipk 
bounded by 2 for k > j. Indeed, in the second case, the density of Dk may be obtained 
by convolution since Dk-i and R,kF{Sk) are independent random variables. 
Therefore, for k < j, P^[\Dk — 0.2| < a\K = j) = l{a>o.2} < 5a since Dfe = for 
k < j and for k > j, P^(|Dfc - 0.2| < a\K = j) = Mv)dv < 4a since Vfe < 2. 

Eventually, 

P^(^fc e W") = P^dA - 0.2| < a) = ^P^(pfc - 0.2| < a\K = j)I>^{K = j) < 5a. 



The results follows. 

Lemma 4.10. For all fc e N, let 

Tjk = {{Tk + u*{Zk)) A Tfc+i) - ((ffc + u*{Zk)) A fk+i) 



□ 



one has for all a > 0, 



\\Vk\\p< ||Tfe-fk||p + 2||Tfc+i -ffe+i||p+ {[u*]^ + 



a 



\\Zk-Zk\\p + 10Cu^a^ 



where [u% = i+gy+^xio^ 

Proof. Let a > 0. Let Ua = [0, 0.2 -a] x {0} x [10"'^; 10"^] x {0}. We will prove that 
the function u*{d, 0, p, 0) is Lipschitz continuous on this set. The function u*{d, 0, p, 0) 
satisfies the following equivalent equations 
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d+dm{u*)^ 0.2 ^ d + p (u* + 1] (^e - ij'j 



= 0.2 



The implicit equation satisfied by u* yields that, on the set Ua, one has u* > 



lO^a. This lower bound will be crucial to prove the Lipschitz continuity. Let d, d' < 
0.2 — a and denote u = u* {d,0, p, 0) and u' = u*{d' ,0, p, 0). Notice that d + dm{u) = 
d' + d„i{u') because they are both worth 0.2. Consequently \dm{u) — djn{u') \ = \d' — d\ 
and, noticing that rj < 2 x 10^ one has 

\d — d'\ ^ p u — u' + rj[e~~^ — e~~) 

> p ^1 — e \u — u'\ 

> 10"'^(l-e~^) \u-u'\ 
that proves the Lipschitz continuity of u* w.r.t. d on 

Similarly, let p, p' £ [10~^; 10~^] and denote u = u*{d, 0, p, 0) and u' = u*{d, 0, p', 0). 
Notice that d + p {u + rj ^e^"? ^ l)) = d + p' (^u' + rj ^e^^ ^ "'^)) ^'^'^^^^'^ they are 
both worth 0.2. Subtracting d + p \u' + rj [e^~ — 1)) in both terms yields 



' + rj{e 1 — e i ) ^ \p — p'\ u' + rj (^e >) — 1^ 



A lower bound for the left-hand side term has already been computed earlier while the 
right hand-side is easily bounded by (C„. -t- 4 x lO'^) |p — p'|, since rj < 2 x 10^, so that 
one has 

(C„. + 4x 10^)|p-p'| > lO-^(l-e-t) \u-u'\ 
that proves the Lipschitz continuity of u* w.r.t. p on Ua- Eventually, for all a > 0, the 



function u* is Lipschitz continuous on Ua with Lipschitz constant [u*]a = 



1+C„« -1-4x10^ 
10-''(l-e~t) ' 



Let fc S N, we now intend to bound ||?yA;||p- Define, as in the proof of Lemma 4.9 
the random variable K = iid{k > such that Gk = 0}, K is the jump happening at 
the end of the protection against the corrosion. 

First, notice that, on the event {k < K} (i.e. when the protection from corrosion is 
still active), one has Zk G for some i G {1, 2, 3} and since the projection defining 
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Zi; from Zj, ensures that they are in the same mode, one has e too. Moreover, 

u*{x) = +00 for all X G -E'(i.i) so that 

lhfcl{fc<K}||p = II {Tk+l - 7A;+i)l{fc<;^}||p < llTfc+i - Tk+l\\p. 

Furthermore, if Zk = A where A denotes the cemetery state, then Z^ — projr^ i^k) — 
A too and one has 77^ = so that 

hfci{fe>/f}llp < hfci{fe>K}i{z,#A}llp <l|7fc - Tk\[p + \\Tk+i - ffe+iiip 

+ ll(w*(^fe) - U*{Zk))l{k>K}'^{Zt^A}\\p- 

Eventually, we intend to bound the last term of the previous sum and we consider 
therefore the event {k > K} n {Zk ^ A}. On the one hand, the random variables Zk 
and Zk both belong to -E(i,o) for some i e {1, 2, 3}. On the other hand, although Um — 
] - 00; 0.2] X for all m e M, one has actually Zk E [0; 0.2] x {0} x [10"'^; lO"""^] x M+ 



p.s. and, according to remark |3.2[ Zk e [0; 0.2] x {0} x [10^^; lO"""^] x M+ p.s. too. 
Combining the two previous remark, one has Zk £ U and Zk £ U where U — [0; 0.2] x 
{0} X [10"'^; 10^5] X {0}. Finally, let a > and notice that C/ C (7^ U [/". One has 

\\{u*{Zk) - u*{Zk))l{k>K}^{z,^A}\\p <A + B 

where 

A=\\{u*{Zk)-u*{Zu))l^^^^-Q^^l{k>K}\\p. 
B ^ Wiu'iZk) - u* {Zk))l{z^eu-'}'^{k>K}\\p- 



The term B is easily bounded thanks to Lemma 4.9 B < 2Cu'P f_i{Zk G U")p < 
lOCu'Ctp. We now turn to the term A and use the Lipschitz continuity of u* on Up 
for any /3 > 0. One has 

A<\\{u*{Zk) -u*{Zk))l^^^^^^^l^^^^^^^l{k>K}\\p 
+ \\{u*{Zk)-u*{Zk))l^^^^jj^^l^-^^^^^l{k>K}\\p 



<[u*].\\Zk-Zk\\p + 2C^4l^^^^~^^l^^ 
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Notice now that 1,^ ~ ^1,^ ~ , < 1,,^ ^ so that finally 

{Z^eUa} {Zk<^Ua} - {|Zfc-Zfc|>f} 



A< [u*]^\\Zk-Zk\\p + 2Cu* (p^{\Zk-Zk\ > 
< [u*]^\\Zk-Zk\\p + ACu' 



\Zk — ZkWp 



and the result follows. 



□ 



The mean exit time Simulation results for the approximation of the mean exit time 
are given m Tabled In order to have a value of reference, a Monte Carlo method (10® 
simulations) yields the value E[T]Monte-Cario — 526 x 10^ h. For the first moment, 
the empirical convergence rate is presented on Figure |5] It is estimated through a 
regression model as —0.38. Remark that it is roughly the same order as the rate of 



convergence of the optimal quantizer (see Theorem 3.11 as here the dimension is 4 



Points in the quantization grids 


PN,1 (Xl03 h) 


relative error to 526 x lO'^h 


20 points 


572 


8.7% 


50 points 


569 


8.2% 


100 points 


557 


5.9% 


200 points 


551 


4.8% 


500 points 


539 


2.5% 



Table 5: Simulations results for the mean exit time 



The exit time distribution Considering the approximation scheme for the exit time 
distribution, one may notice that the quantized value PAr(s) is not necessary smaller 
than 1. Therefore, it appears natural to replace PAr(s) by pAf(s) A 1. This does not 
change the convergence theorem and can only improve the approximation error. It is 
equally possible, and this is done in the results below, to replace Pn{s) by Z"^^^ since 

Pn(0) 

Pn{0) goes to 1. 

Figure [6] presents the survival function of r obtained through Monte Carlo Simula- 



Exit time for PDMP 



41 




Figure 5: Logarithm of the error w.r.t. the logarithm of the number of points in the 
quantization grids for the first moment of the corrosion process. 



tions (the dashed line), through our approximation scheme (the solid line) and the error. 
Table [6] contains the empirical error for different numbers of points in the quantization 
grids. For the survival function, the empirical convergence rate is presented on Figure 
[7] It is estimated through a regression model as —0.63. Remark that it is roughly the 



same order as the rate of convergence of the optimal quantizer (see Theorem 3.11 as 
here the dimension is 4. 



Points in the quantization grids 


maxs \pn{s) -pn{s)\ 


20 points 


0.145 


50 points 


0.119 


100 points 


0.040 


200 points 


0.039 


500 points 


0.020 



Table 6: Simulations results for the distribution 



The convergence of the approximation scheme in the corrosion model appears to be 
slightly slower than in the previous example. This is due to the higher dimension of 
the process that has to be quantized, which is 4 in the case of the corrosion model and 
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Figure 6: Survival function of r obtained through Monte Carlo simulations (dashed), 
quantized approximation (solid) and the error with 500 points in the quantization grids. 

1 in the case of the Poisson process. 

Remark 4.11. By using Monte Carlo simulations, we can approximate the value of q. 
One has q ~ 0.0187 for 10^ histories. 



5. Advantages and practical interest of our approach 

Let us describe the practical interest of our approach. 

• The quantizations grids only have to be computed once and for all and can be 
used for several purposes. Moreover, once they are obtained, the procedures 
leading to pn{s) and to pnj can be achieved very simply since we only have to 
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Figure 7: Logarithm of the error w.r.t. the logarithm of the number of points in the 
quantization grids for the survival of the corrosion process. 



compute finite sums. 

• Concerning the distribution, since pn{s) can be computed almost instantly for 
any value of s, the whole survival function can be obtained very quickly. Similarly, 
concerning the moments, PN,j can be computed very quickly for any j, so that 
any moment is almost instantly available. 

• Furthermore, in both cases, one may decide to change the set U and consider 
the exit time r' from a new set U' . This will yield new sequences {qk)k, (j^k.j)k 
and {pk,j)k in the case of the j-th moment approximation or new sequences {qk)k, 
{rk{s))k and {pk{s))k if we are interested in the distribution. These new sequences 
are obtained quickly and easily since the quantized process remains the same and 
we only have to compute finite sums. Of course, the set U' must be such that 



Assumptions 2.3 to 2.7 remain true and such that P^t(T;v < t') remains small 
without changing the computation horizon N. This last condition is fulfilled if, 
for instance, [/' C U. This flexibility is an important advantage of our method 
over, for instance, a Monte Carlo method. 
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